setwd(dirname(rstudioapi::getActiveDocumentContext()$path))

## Import packages
require(cregg)
require(ggplot2)

#################
## Import Data ##
#################

## Experiment 1
dcj_exp1 <- readRDS("elecvspref_jan22_conjoint_v4.1.rds")

## Experiment 2
dcj_exp2 <- readRDS("elecvspref_mar16_conjoint_v4.1.rds")

##################
## Experiment 1 ##
##################

## Clean Data
dim(dcj_exp1)
## Experimental Assignment Mechanically Missing
dcj_exp1 <- subset(dcj_exp1, !is.na(dcj_exp1$candgender))
dim(dcj_exp1)
## Non-Responses
dcj_exp1 <- subset(dcj_exp1, !is.na(dcj_exp1$selected))
dim(dcj_exp1)
(21624-21218)/21624 # Nonresponse rate

## Adjust Variables
dcj_exp1$dv <- 
  factor(ifelse(dcj_exp1$electability==1, "Expectation", "Preference"),
         levels = c("Preference", "Expectation"))
dcj_exp1$candgender <- 
  factor(gsub("^\\(.*\\) ","", dcj_exp1$candgender),
         levels = gsub("^\\(.*\\) ","", levels(dcj_exp1$candgender)))
dcj_exp1$candparty <- 
  factor(gsub("^\\(.*\\) ","", dcj_exp1$candparty),
         levels = gsub("^\\(.*\\) ","", levels(dcj_exp1$candparty)))
dcj_exp1$candage <- 
  factor(gsub("^\\(.*\\) ","", dcj_exp1$candage),
         levels = gsub("^\\(.*\\) ","", levels(dcj_exp1$candage)))
dcj_exp1$candexperience <- 
  factor(gsub("^\\(.*\\) ","", dcj_exp1$candexperience),
         levels = gsub("^\\(.*\\) ","", levels(dcj_exp1$candexperience)))
dcj_exp1$candedu <- 
  factor(gsub("^\\(.*\\) ","", dcj_exp1$candedu),
         levels = gsub("^\\(.*\\) ","", levels(dcj_exp1$candedu)))
dcj_exp1$candmarriage <- 
  factor(gsub("^\\(.*\\) ","", dcj_exp1$candmarriage),
         levels = gsub("^\\(.*\\) ","", levels(dcj_exp1$candmarriage)))
dcj_exp1$candkidage <- 
  factor(gsub("^\\(.*\\) ","", dcj_exp1$candkidage),
         levels = gsub("^\\(.*\\) ","", levels(dcj_exp1$candkidage)))
dcj_exp1$candlivetype <- 
  factor(gsub("^\\(.*\\) ","", dcj_exp1$candlivetype),
         levels = gsub("^\\(.*\\) ","", levels(dcj_exp1$candlivetype)))

## Variable Labels 
vlab_exp1 <- list("candgender"="Gender",
             "candparty"="Party Affliation",
             "candage"="Age",
             "candexperience"="Political Experience",
             "candedu"="Educational Attainment",
             "candmarriage"="Marital Status",
             "candkidage"="Youngest Child",
             "candlivetype"="Residence")

## Marginal Means ##

resexp1_each <- cj(selected ~ candgender + candparty + candage + 
                 candexperience + candedu + candmarriage + 
                 candkidage + candlivetype, 
               data = dcj_exp1,
               feature_labels = vlab_exp1,
               by = ~dv, 
               id = ~Response.ID, estimate="mm", h0=0.5)
head(resexp1_each)
resexp1_diff <- cj(selected ~ candgender + candparty + candage + 
                 candexperience + candedu + candmarriage + 
                 candkidage + candlivetype, 
               data = dcj_exp1,
               feature_labels = vlab_exp1,
               by = ~dv, 
               id = ~Response.ID, estimate="mm_diff")
head(resexp1_diff)

resexp1_both <- rbind(resexp1_each,resexp1_diff)
resexp1_both$BY <- 
  factor(resexp1_both$BY, 
         levels = c("Preference","Expectation",              
                    "Expectation - Preference"))
resexp1_both$level <- factor(resexp1_both$level, 
                         levels = rev(levels(resexp1_both$level)))

setfacetexp1_vline <- 
  data.frame(BY = factor(rep(c("Preference","Expectation",              
                               "Expectation - Preference"), 
                             each = 2),
                         levels = 
                           c("Preference","Expectation",              
                             "Expectation - Preference")),
             VL = rep(c(0.5,0.5,0),each=2))
setfacetexp1_lims <- 
  data.frame(BY = factor(rep(c("Preference","Expectation",              
                               "Expectation - Preference"), 
                             each = 2),
                         levels = 
                           c("Preference","Expectation",              
                             "Expectation - Preference")),
             level = factor("Male", levels(resexp1_both$level)), 
             feature = factor("Gender", levels(resexp1_both$feature)),
             lims = c(0.3,0.7,0.3,0.7,-0.15,0.15))

## Plot!
## Figure A1 (Appendix B)
require(ggplot2)
p_resexp1 <- 
  ggplot(resexp1_both) + 
  geom_vline(data=setfacetexp1_vline, aes(xintercept=VL), 
             linetype=2, color="gray70") + 
  geom_errorbarh(aes(xmin=lower,xmax=upper,y=level), height=0) +
  geom_point(aes(x=estimate,y=level), size=1) + 
  geom_point(data=setfacetexp1_lims, aes(x=lims,y=level), 
             size=1, alpha=0) + 
  facet_grid(feature~BY, scales = "free", space = "free") + 
  labs(x="Marginal Means Estimates (with 95% Confidence Intervals)",
       y=NULL) + 
  theme_bw(base_family = "serif") + 
  theme(strip.text.y = element_text(angle=0),
        plot.title = element_text(hjust=0.5), 
        legend.position = "bottom")
p_resexp1
ggsave("resexp1_v4.pdf", width=8, height=7)

### Extracting p-values
### Figure A2 (Appendix B)
require(ggplot2)
p_resexp1_pval <- 
  ggplot(resexp1_both) + 
  geom_text(aes(x=0.5,
                y=level,label=sprintf("%.03f",p)), size=3) + 
  geom_point(aes(x=0.48,y=level), size=1, alpha=0) + 
  geom_point(aes(x=0.52,y=level), size=1, alpha=0) + 
  facet_grid(feature~BY, scales = "free", space = "free") + 
  labs(x="p-values from significance test",
       y=NULL,
       caption="Note: H0=0.5 for preference and expectation tasks. H0=0 for preference-expecation gap.") + 
  theme_minimal(base_family = "serif") + 
  theme(strip.text.y = element_text(angle=0),
        plot.title = element_text(hjust=0.5), 
        legend.position = "bottom",
        axis.text.x = element_blank())
p_resexp1_pval
ggsave("resexp1_pval_v4.pdf", width=8, height=7)

## AMCE (only estimation) ##

resexp1_amce_each <- cj(selected ~ candgender + candparty + candage + 
                     candexperience + candedu + candmarriage + 
                     candkidage + candlivetype, 
                   data = dcj_exp1,
                   feature_labels = vlab_exp1,
                   by = ~dv, 
                   id = ~Response.ID, estimate="amce")
head(resexp1_amce_each)
resexp1_amce_diff <- cj(selected ~ candgender + candparty + candage + 
                     candexperience + candedu + candmarriage + 
                     candkidage + candlivetype, 
                   data = dcj_exp1,
                   feature_labels = vlab_exp1,
                   by = ~dv, 
                   id = ~Response.ID, estimate="amce_diff")
head(resexp1_amce_diff)

resexp1_amce_both <- rbind(resexp1_amce_each,resexp1_amce_diff)
resexp1_amce_both$BY <- 
  factor(resexp1_amce_both$BY, 
         levels = c("Preference","Expectation",              
                    "Expectation - Preference"))
resexp1_amce_both$level <- factor(resexp1_amce_both$level, 
                             levels = rev(levels(resexp1_amce_both$level)))

##################
## Experiment 2 ##
##################

## Clean Data
dim(dcj_exp2)
## Experimental Assignment Mechanically Missing
dcj_exp2 <- subset(dcj_exp2, !is.na(dcj_exp2$candgender))
dim(dcj_exp2)
## Nonresponses
dcj_exp2 <- subset(dcj_exp2, !is.na(dcj_exp2$selected))
dim(dcj_exp2)
(38464-38222)/38464 # Nonresponse rate

## Adjust Variables
dcj_exp2$dv <- 
  factor(ifelse(dcj_exp2$electability==1, "Expectation", "Preference"),
         levels = c("Preference", "Expectation"))
dcj_exp2$candgender <- 
  factor(gsub("^\\(.*\\) ","", dcj_exp2$candgender),
         levels = gsub("^\\(.*\\) ","", levels(dcj_exp2$candgender)))
dcj_exp2$candparty <- 
  factor(gsub("^\\(.*\\) ","", dcj_exp2$candparty),
         levels = gsub("^\\(.*\\) ","", levels(dcj_exp2$candparty)))
dcj_exp2$candage <- 
  factor(gsub("^\\(.*\\) ","", dcj_exp2$candage),
         levels = gsub("^\\(.*\\) ","", levels(dcj_exp2$candage)))
dcj_exp2$candexperience <- 
  factor(gsub("^\\(.*\\) ","", dcj_exp2$candexperience),
         levels = gsub("^\\(.*\\) ","", levels(dcj_exp2$candexperience)))
dcj_exp2$candedu <- 
  factor(gsub("^\\(.*\\) ","", dcj_exp2$candedu),
         levels = gsub("^\\(.*\\) ","", levels(dcj_exp2$candedu)))
dcj_exp2$candfam <- 
  factor(gsub("^\\(.*\\) ","", dcj_exp2$candfam),
         levels = gsub("^\\(.*\\) ","", levels(dcj_exp2$candfam)))
dcj_exp2$candlivetype <- 
  factor(gsub("^\\(.*\\) ","", dcj_exp2$candlivetype),
         levels = gsub("^\\(.*\\) ","", levels(dcj_exp2$candlivetype)))
dcj_exp2$candorigin <- 
  factor(gsub("^\\(.*\\) ","", dcj_exp2$candorigin),
         levels = gsub("^\\(.*\\) ","", levels(dcj_exp2$candorigin)))
dcj_exp2$candimppol <- 
  factor(gsub("^\\(.*\\) ","", dcj_exp2$candimppol),
         levels = gsub("^\\(.*\\) ","", levels(dcj_exp2$candimppol)))

## Variable Labels 
vlab_exp2 <- list("candgender"="Gender",
                  "candparty"="Party Affliation",
                  "candage"="Age",
                  "candexperience"="Political Experience",
                  "candedu"="Educational Attainment",
                  "candfam"="Family",
                  "candorigin"="Native Place",
                  "candlivetype"="Residence",
                  "candimppol"="Policy Focus")

## house (Marginal Means) ##

resexp2_each_house <- cj(selected ~ candgender + candparty + candage + 
                           candexperience + candedu + candfam + 
                           candlivetype + candorigin + candimppol, 
                         data = subset(dcj_exp2, conjoint_eleclevel=="house"),
                         feature_labels = vlab_exp2,
                         by = ~dv, 
                         id = ~Response.ID, estimate="mm", h0=0.5)
head(resexp2_each_house)
resexp2_diff_house <- cj(selected ~ candgender + candparty + candage + 
                           candexperience + candedu + candfam + 
                           candlivetype + candorigin + candimppol, 
                         data = subset(dcj_exp2, conjoint_eleclevel=="house"),
                         feature_labels = vlab_exp2,
                         by = ~dv, 
                         id = ~Response.ID, estimate="mm_diff")
head(resexp2_diff_house)

resexp2_both_house <- rbind(resexp2_each_house,
                            resexp2_diff_house)
resexp2_both_house$BY <- 
  factor(resexp2_both_house$BY, 
         levels = c("Preference","Expectation",              
                    "Expectation - Preference"))
resexp2_both_house$level <- 
  factor(resexp2_both_house$level, 
         levels = rev(levels(resexp2_both_house$level)))

setfacetexp2_vline_house <- 
  data.frame(BY = factor(rep(c("Preference","Expectation",              
                               "Expectation - Preference"), 
                             each = 2),
                         levels = 
                           c("Preference","Expectation",              
                             "Expectation - Preference")),
             VL = rep(c(0.5,0.5,0),each=2))
setfacetexp2_lims_house <- 
  data.frame(BY = factor(rep(c("Preference","Expectation",              
                               "Expectation - Preference"), 
                             each = 2),
                         levels = 
                           c("Preference","Expectation",              
                             "Expectation - Preference")),
             level = factor("Male", levels(resexp2_both_house$level)), 
             feature = factor("Gender", levels(resexp2_both_house$feature)),
             lims = c(0.3,0.7,0.3,0.7,-0.15,0.15))

## Plot!
## Figure A3 (Appendix B)
require(ggplot2)
p_resexp2_house <- 
  ggplot(resexp2_both_house) + 
  geom_vline(data=setfacetexp2_vline_house, aes(xintercept=VL), 
             linetype=2, color="gray70") + 
  geom_errorbarh(aes(xmin=lower,xmax=upper,y=level), height=0) +
  geom_point(aes(x=estimate,y=level), size=1) + 
  geom_point(data=setfacetexp2_lims_house, aes(x=lims,y=level), 
             size=1, alpha=0) + 
  facet_grid(feature~BY, scales = "free", space = "free") + 
  labs(x="Marginal Means Estimates (with 95% Confidence Intervals)",
       y=NULL) + 
  theme_bw(base_family = "serif") + 
  theme(strip.text.y = element_text(angle=0),
        plot.title = element_text(hjust=0.5), 
        legend.position = "bottom")
p_resexp2_house
ggsave("resexp2_house_v4.pdf", width=8, height=7)

### Extracting p-values
### Figure A4 (Appendix B)
require(ggplot2)
p_resexp2_house_pval <- 
  ggplot(resexp2_both_house) + 
  geom_text(aes(x=0.5,
                y=level,label=sprintf("%.03f",p)), size=3) + 
  geom_point(aes(x=0.48,y=level), size=1, alpha=0) + 
  geom_point(aes(x=0.52,y=level), size=1, alpha=0) + 
  facet_grid(feature~BY, scales = "free", space = "free") + 
  labs(x="p-values from significance test",
       y=NULL,
       caption="Note: H0=0.5 for preference and expectation tasks. H0=0 for preference-expecation gap.") + 
  theme_minimal(base_family = "serif") + 
  theme(strip.text.y = element_text(angle=0),
        plot.title = element_text(hjust=0.5), 
        legend.position = "bottom",
        axis.text.x = element_blank())
p_resexp2_house_pval
ggsave("resexp2_house_pval_v4.pdf", width=8, height=7)

## house (AMCE, only estimation) ##

resexp2_amce_each_house <- cj(selected ~ candgender + candparty + candage + 
                           candexperience + candedu + candfam + 
                           candlivetype + candorigin + candimppol, 
                         data = subset(dcj_exp2, conjoint_eleclevel=="house"),
                         feature_labels = vlab_exp2,
                         by = ~dv, 
                         id = ~Response.ID, estimate="amce")
head(resexp2_amce_each_house)
resexp2_amce_diff_house <- cj(selected ~ candgender + candparty + candage + 
                           candexperience + candedu + candfam + 
                           candlivetype + candorigin + candimppol, 
                         data = subset(dcj_exp2, conjoint_eleclevel=="house"),
                         feature_labels = vlab_exp2,
                         by = ~dv, 
                         id = ~Response.ID, estimate="amce_diff")
head(resexp2_amce_diff_house)

resexp2_amce_both_house <- rbind(resexp2_amce_each_house,
                            resexp2_amce_diff_house)
resexp2_amce_both_house$BY <- 
  factor(resexp2_amce_both_house$BY, 
         levels = c("Preference","Expectation",              
                    "Expectation - Preference"))
resexp2_amce_both_house$level <- 
  factor(resexp2_amce_both_house$level, 
         levels = rev(levels(resexp2_amce_both_house$level)))

## municipal (Marignal Means) ##

## Calculate Marginal Means for each minister
resexp2_each_municipal <- cj(selected ~ candgender + candparty + candage + 
                           candexperience + candedu + candfam + 
                           candlivetype + candorigin + candimppol, 
                         data = subset(dcj_exp2, conjoint_eleclevel=="municipal"),
                         feature_labels = vlab_exp2,
                         by = ~dv, 
                         id = ~Response.ID, estimate="mm", h0=0.5)
head(resexp2_each_municipal)
resexp2_diff_municipal <- cj(selected ~ candgender + candparty + candage + 
                           candexperience + candedu + candfam + 
                           candlivetype + candorigin + candimppol, 
                         data = subset(dcj_exp2, conjoint_eleclevel=="municipal"),
                         feature_labels = vlab_exp2,
                         by = ~dv, 
                         id = ~Response.ID, estimate="mm_diff")
head(resexp2_diff_municipal)

resexp2_both_municipal <- rbind(resexp2_each_municipal,
                            resexp2_diff_municipal)
resexp2_both_municipal$BY <- 
  factor(resexp2_both_municipal$BY, 
         levels = c("Preference","Expectation",              
                    "Expectation - Preference"))
resexp2_both_municipal$level <- 
  factor(resexp2_both_municipal$level, 
         levels = rev(levels(resexp2_both_municipal$level)))

setfacetexp2_vline_municipal <- 
  data.frame(BY = factor(rep(c("Preference","Expectation",              
                               "Expectation - Preference"), 
                             each = 2),
                         levels = 
                           c("Preference","Expectation",              
                             "Expectation - Preference")),
             VL = rep(c(0.5,0.5,0),each=2))

setfacetexp2_lims_municipal <- 
  data.frame(BY = factor(rep(c("Preference","Expectation",              
                               "Expectation - Preference"), 
                             each = 2),
                         levels = 
                           c("Preference","Expectation",              
                             "Expectation - Preference")),
             level = factor("Male", levels(resexp2_both_municipal$level)), 
             feature = factor("Gender", levels(resexp2_both_municipal$feature)),
             lims = c(0.3,0.7,0.3,0.7,-0.15,0.15))

## Plot!
## Figure A5 (Appendix B)
require(ggplot2)
p_resexp2_municipal <- 
  ggplot(resexp2_both_municipal) + 
  geom_vline(data=setfacetexp2_vline_municipal, aes(xintercept=VL), 
             linetype=2, color="gray70") + 
  geom_errorbarh(aes(xmin=lower,xmax=upper,y=level), height=0) +
  geom_point(aes(x=estimate,y=level), size=1) + 
  geom_point(data=setfacetexp2_lims_municipal, aes(x=lims,y=level), 
             size=1, alpha=0) + 
  facet_grid(feature~BY, scales = "free", space = "free") + 
  labs(x="Marginal Means Estimates (with 95% Confidence Intervals)",
       y=NULL) + 
  theme_bw(base_family = "serif") + 
  theme(strip.text.y = element_text(angle=0),
        plot.title = element_text(hjust=0.5), 
        legend.position = "bottom")
p_resexp2_municipal
ggsave("resexp2_municipal_v4.pdf", width=8, height=7)

### Extracting p-values
### Figure A6 (Appendix B)
require(ggplot2)
p_resexp2_municipal_pval <- 
  ggplot(resexp2_both_municipal) + 
  geom_text(aes(x=0.5,
                y=level,label=sprintf("%.03f",p)), size=3) + 
  geom_point(aes(x=0.48,y=level), size=1, alpha=0) + 
  geom_point(aes(x=0.52,y=level), size=1, alpha=0) + 
  facet_grid(feature~BY, scales = "free", space = "free") + 
  labs(x="p-values from significance test",
       y=NULL,
       caption="Note: H0=0.5 for preference and expectation tasks. H0=0 for preference-expecation gap.") + 
  theme_minimal(base_family = "serif") + 
  theme(strip.text.y = element_text(angle=0),
        plot.title = element_text(hjust=0.5), 
        legend.position = "bottom",
        axis.text.x = element_blank())
p_resexp2_municipal_pval
ggsave("resexp2_municipal_pval_v4.pdf", width=8, height=7)

## municipal (AMCE, only estimation) ##

## Calculate Marginal Means for each minister
resexp2_amce_each_municipal <- cj(selected ~ candgender + candparty + candage + 
                               candexperience + candedu + candfam + 
                               candlivetype + candorigin + candimppol, 
                             data = subset(dcj_exp2, conjoint_eleclevel=="municipal"),
                             feature_labels = vlab_exp2,
                             by = ~dv, 
                             id = ~Response.ID, estimate="amce")
head(resexp2_amce_each_municipal)
resexp2_amce_diff_municipal <- cj(selected ~ candgender + candparty + candage + 
                               candexperience + candedu + candfam + 
                               candlivetype + candorigin + candimppol, 
                             data = subset(dcj_exp2, conjoint_eleclevel=="municipal"),
                             feature_labels = vlab_exp2,
                             by = ~dv, 
                             id = ~Response.ID, estimate="amce_diff")
head(resexp2_amce_diff_municipal)

resexp2_amce_both_municipal <- rbind(resexp2_amce_each_municipal,
                                resexp2_amce_diff_municipal)
resexp2_amce_both_municipal$BY <- 
  factor(resexp2_amce_both_municipal$BY, 
         levels = c("Preference","Expectation",              
                    "Expectation - Preference"))
resexp2_amce_both_municipal$level <- 
  factor(resexp2_amce_both_municipal$level, 
         levels = rev(levels(resexp2_amce_both_municipal$level)))

##############
## Figure 1 ##
##############

tmp1 <- subset(resexp1_both, feature=="Gender")
tmp1$exp <- "Experiment 1\n[Without Policy Attribute]\n(House of Representatives)"
tmp2a <- subset(resexp2_both_house, feature=="Gender")
tmp2a$exp <- "Experiment 2\n[With Policy Attribute]\n(House of Representatives)"
tmp2b <- subset(resexp2_both_municipal, feature=="Gender")
tmp2b$exp <- "Experiment 2\n[With Policy Attribute]\n(Municipal Council)"

fig1dt <- rbind(tmp1, tmp2a, tmp2b)
fig1dt$exp <- factor(fig1dt$exp, unique(fig1dt$exp))

setfacetfig1_vline <- 
  data.frame(BY = factor(rep(c("Preference","Expectation",              
                               "Expectation - Preference"), 
                             each = 2),
                         levels = 
                           c("Preference","Expectation",              
                             "Expectation - Preference")),
             VL = rep(c(0.5,0.5,0),each=2))
setfacetfig1_lims <- 
  data.frame(BY = factor(rep(c("Preference","Expectation",              
                               "Expectation - Preference"), 
                             each = 2),
                         levels = 
                           c("Preference","Expectation",              
                             "Expectation - Preference")),
             level = factor("Male", levels(fig1dt$level)), 
             feature = factor("Gender", levels(fig1dt$feature)),
             lims = c(0.45,0.55,0.45,0.55,-0.06,0.06))

## Plot!
require(ggplot2)
p_resfig1 <- 
  ggplot(fig1dt) + 
  geom_vline(data=setfacetfig1_vline, aes(xintercept=VL), 
             linetype=2, color="gray70") + 
  geom_errorbarh(aes(xmin=lower,xmax=upper,y=level), height=0) +
  geom_point(aes(x=estimate,y=level), size=1.5) + 
  geom_point(data=setfacetfig1_lims, aes(x=lims,y=level), 
             size=1, alpha=0) + 
  facet_grid(exp~BY, scales = "free", 
             switch = "y") + 
  scale_x_continuous(breaks = c(-0.05,0,0.05,0.47,0.5,0.53)) + 
  labs(x="Marginal Means Estimates (with 95% Confidence Intervals)",
       y=NULL) + 
  theme_bw(base_family = "serif") + 
  theme(plot.title = element_text(hjust=0.5), 
        legend.position = "bottom",
        panel.grid = element_blank(),
        strip.placement = "outside",
        strip.text.y.left = element_text(angle=0),
        strip.background.y = element_blank())
p_resfig1
ggsave("fig1_v4.pdf", width=7, height=3.5)

####################################
## Appendix C: Figure 1 AMCE ver. ##
####################################

tmp1 <- subset(resexp1_amce_both, feature=="Gender")
tmp1$exp <- "Experiment 1\n[Without Policy Attribute]\n(House of Representatives)"
tmp2a <- subset(resexp2_amce_both_house, feature=="Gender")
tmp2a$exp <- "Experiment 2\n[With Policy Attribute]\n(House of Representatives)"
tmp2b <- subset(resexp2_amce_both_municipal, feature=="Gender")
tmp2b$exp <- "Experiment 2\n[With Policy Attribute]\n(Municipal Council)"

fig1_amcedt <- rbind(tmp1, tmp2a, tmp2b)
fig1_amcedt$exp <- factor(fig1_amcedt$exp, unique(fig1_amcedt$exp))

setfacetfig1_amce_vline <- 
  data.frame(BY = factor(rep(c("Preference","Expectation",              
                               "Expectation - Preference"), 
                             each = 2),
                         levels = 
                           c("Preference","Expectation",              
                             "Expectation - Preference")),
             VL = rep(c(0,0,0),each=2))
setfacetfig1_amce_lims <- 
  data.frame(BY = factor(rep(c("Preference","Expectation",              
                               "Expectation - Preference"), 
                             each = 2),
                         levels = 
                           c("Preference","Expectation",              
                             "Expectation - Preference")),
             level = factor("Male", levels(fig1_amcedt$level)), 
             feature = factor("Gender", levels(fig1_amcedt$feature)),
             lims = c(-0.08,0.08,-0.08,0.08,-0.14,0.02))
fig1_amcedt

## Plot!
## Figure A7 (Appendix C)
require(ggplot2)
p_resfig1_amce <- 
  ggplot(fig1_amcedt) + 
  geom_vline(data=setfacetfig1_amce_vline, aes(xintercept=VL), 
             linetype=2, color="gray70") + 
  geom_errorbarh(aes(xmin=lower,xmax=upper,y=level), height=0) +
  geom_point(aes(x=estimate,y=level), size=1.5) + 
  geom_point(data=setfacetfig1_amce_lims, aes(x=lims,y=level), 
             size=1, alpha=0) + 
  facet_grid(exp~BY, scales = "free", 
             switch = "y") + 
  scale_x_continuous(breaks = c(-0.1,-0.05,0,0.05)) + 
  labs(x="Average Marginal Component Effect (with 95% Confidence Intervals)",
       y=NULL) + 
  theme_bw(base_family = "serif") + 
  theme(plot.title = element_text(hjust=0.5), 
        legend.position = "bottom",
        panel.grid = element_blank(),
        strip.placement = "outside",
        strip.text.y.left = element_text(angle=0),
        strip.background.y = element_blank())
p_resfig1_amce
ggsave("fig1_amce_v4.pdf", width=7, height=3.5)

### Extracting p-values
### Figure A8 (Appendix C)
require(ggplot2)
p_resfig1_amce_pval <- 
  ggplot(fig1_amcedt) + 
  geom_text(aes(x=0.5,
                y=level,label=sprintf("%.03f",p)), size=3) + 
  geom_point(aes(x=0.48,y=level), size=1, alpha=0) + 
  geom_point(aes(x=0.52,y=level), size=1, alpha=0) + 
  facet_grid(exp~BY, scales = "free", 
             switch = "y") + 
  labs(x="p-values from significance test (H0=0)",
       y=NULL) + 
  theme_minimal(base_family = "serif") + 
  theme(plot.title = element_text(hjust=0.5), 
        legend.position = "bottom",
        strip.placement = "outside",
        strip.text.y.left = element_text(angle=0),
        strip.background.y = element_blank(),
        axis.text.x = element_blank())
p_resfig1_amce_pval
ggsave("fig1_amce_pval_v4.pdf", width=7, height=3.5)
